Nature of the anomalies in the supercooled liquid state of the 
mW model of water 

Vincent Holten, 1 David T. Limmer, 2 Valeria Molinero, 3 and Mikhail A. Anisimov 1 ' a) 

1 ^Institute for Physical Science and Technology and Department of Chemical and Biomolecular Engineering, 
University of Maryland, College Park, Maryland 20742, USA 

^Department of Chemistry, University of California, Berkeley, California 94720, USA 
^Department of Chemistry, University of Utah, Salt Lake City, Utah 84112-0580, USA 

(Dated: 3 April 2013) 

The thermodynamic properties of the supercooled liquid state of the mW model of water show anoma- 
lous behavior. Like in real water, the heat capacity and compressibility sharply increase upon super- 
cooling. One of the possible explanations of these anomalies, the existence of a second (liquid-liquid) 
critical point, is not supported by simulations for this model. In this work, we reproduce the anomalies 
of the mW model with two thermodynamic scenarios: one based on a non-ideal "mixture" with two 
different types of local order of the water molecules, and one based on weak crystallization theory. 
We show that both descriptions accurately reproduce the model's basic thermodynamic properties. 
However, the coupling constant required for the power laws implied by weak crystallization theory 
is too large relative to the regular backgrounds, contradicting assumptions of weak crystallization 
theory. Fluctuation corrections outside the scope of this work would be necessary to fit the forms 
predicted by weak crystallization theory. For the two-state approach, the direct computation of the 
low-density fraction of molecules in the mW model is in agreement with the prediction of the phe- 
nomenological equation of state. The non-ideality of the "mixture" of the two states never becomes 
strong enough to cause liquid-liquid phase separation, also in agreement with simulation results. 



I. INTRODUCTION 

The peculiar properties of supercooled water continue to 
gain interest. In the supercooled region, the thermody- 
namic response functions, namely heat capacity, 1 thermal 
expansivity, 2,3 and compressibility, 4 show strong temperature 
dependences suggesting a possible divergence at a tempera- 
ture just below the homogeneous ice nucleation limit. One of 
the scenarios to explain the anomalous behavior of real water 
is the existence of a liquid-liquid transition terminated by the 
liquid-liquid critical point. 5 " 14 An explicit equation of state 
based on this scenario is able to accurately represent all exper- 
imental data on the thermodynamic properties of supercooled 
water. 15 If it exists, the second critical point of water cannot 
be directly observed in a bulk experiment, because it is located 
in "no man's land," below the homogeneous ice nucleation 
temperature. 16 ' 17 Computer simulations of water can provide 
additional insights into the nature of water's anomalies. The 
mW model devised by Molinero and Moore 18 represents the 
water molecule as a single atom with only short-range in- 
teractions, and is is suitable for fast computations. The mW 
model imitates the anomalous behavior of cold and super- 
cooled water, including the density maximum and the increase 
of the heat capacity and compressibility in the supercooled 
region. 18 " 20 Using molecular dynamics simulations, Limmer 
and Chandler 20 have shown that the mW model does not ex- 
hibit a second critical point or liquid-liquid separation in the 



range studied (0 MPa to 290 MPa, down to 170 K). Indeed, 
Moore and Molinero 19 have demonstrated that in this model 
the supercooled liquid can no longer be equilibrated before it 
crystallizes and there is no sign of a liquid-liquid transition 
at supercooling. This raises the question: what is the origin of 
the anomalies in the mW model? In this work, we explain and 
reproduce these anomalies with a thermodynamic equation of 
state based on a non-ideal "mixture" of two kinds of molecu- 
lar environments, in which the non-ideality is mainly entropy 
driven. By analyzing experimental data with this equation 
of state, Anisimov and coworkers 15 ' 21 have concluded that a 
liquid-liquid transition can occur if effects of crystallization 
are neglected. We show that for the mW model the nonideal- 
ity of the free energy of mixing never becomes large enough to 
cause liquid-liquid phase separation. Finally, we show that the 
properties of the mW model may be described by power laws 
suggested by weak crystallization theory, which predicts ap- 
parently diverging corrections to the regular thermodynamic 
properties as a result of fluctuations in the translational order 
parameter, close to the limit of stability of the liquid phase. 
However, the resulting value of the coupling constant is large, 
which contradicts the basic assumption of the theory that the 
corrections must always remain smaller than the regular back- 
grounds, suggesting that fluctuation corrections beyond what 
are considered here may be important. 



II. TWO-STATE THERMODYNAMICS OF LIQUID 
WATER 
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We assume liquid water at low temperatures to be a mixture of 
two interconvertible states or structures, a high-density state A 



2 



and a low-density state B. The fraction of molecules in state 
B is denoted by x, and is controlled by the 'reaction' 



B. 



(1) 



The states A and B could correspond to different arrangements 
of the hydrogen-bonded network. 22 In water-like atomistic 
models, such as the mW model, these states could correspond 
to two kinds of local coordination of the water molecules. 23 

In a real liquid, molecular configurations form a continu- 
ous distribution of coordination numbers and local structures. 
Therefore, for real water the division of molecular configura- 
tions into two states seems a gross simplification. However, 
such an approach could serve as a first approximation if the 
distribution can be decomposed in two populations with dis- 
tinct properties. A similar concept, "quasi-binary approxima- 
tion," is commonly used to describe the properties of multi- 
component fluids. 24 As for any phenomenological model, the 
utility of the two-state approximation is to be provided by a 
comparison with experimental or computational data. 

Two-state equations of state have become popular to ex- 
plain liquid polyamorphism. 25-28 Ponyatovsky et al. 29 and 
Moynihan 30 assumed that water could be considered as a 'reg- 
ular binary solution' of two states, which implies that the 
phase separation is driven by energy, and they qualitatively 
reproduced the thermodynamic anomalies of water. Cuthbert- 
son and Poole 31 showed that the ST2 water model can be de- 
scribed by a regular-solution two-state equation. Bertrand and 
Anisimov 21 introduced a two-state equation of state where 
water is assumed to be an athermal solution, which undergoes 
phase separation driven by non-ideal entropy upon increase of 
the pressure. This equation of state was used by Holten and 
Anisimov 15 to successfully describe the properties of real su- 
percooled water. 

In general, the molar Gibbs energy of a two-state mixture is 

G = G A +xG BA +RT[x\nx+ (1 -x)ln(l -x)] + G E , (2) 

where R is the gas constant, T is the temperature, G = 
qB qA j s tne (}iff ererice m molar Gibbs energy between pure 
configurations A and B, and G E is the excess Gibbs energy of 
mixing. The difference G BA is related to the equilibrium con- 
stant K of "reaction" (1) as 



\nK(T,P) 



-BA 



RT 



(3) 



where P is the pressure. For the application to the mW model, 
we adopt a linear expression for InK as the simplest approxi- 
mation. 



with 



lnK = A (AT + aP), (4) 



Af = (T - T ) /T , P = P/poRTo, (5) 



where 7b is the temperature at which \nK = for zero pres- 
sure, and po is a reference density. The parameter a is pro- 
portional to the slope of the InK = line in the T—P phase 



diagram. The parameter A is proportional to the heat of reac- 
tion (1), while the product v = Xa is proportional to the vol- 
ume change of the reaction. 
The excess Gibbs energy, 



G E =H E -TS E , 



(6) 



causes the non-ideality of the mixture and is the sum of con- 
tributions of the enthalpy of mixing H E and excess entropy 
S E . 

The Gibbs energy G A of the pure structure A defines the 
background of the properties and is approximated as 



G A =RToY J c m „^f) m P", 



(7) 



where m and n are integers and c m „ are adjustable coefficients. 



Regular solution 

In the case of a regular solution, the non-ideality is entirely 
associated with the enthalpy of mixing. Taking H E = wx{ \ — 
x), we obtain 



G = G A + xG aA + RT [xlnx + ( 1 - x) ln( 1 - x)] + wx{ 1 



-*). 
(8) 

The interaction parameter w is assumed to be independent 
of temperature, but may depend on pressure. Considering 
x, the fraction of B, as the reaction coordinate or extent of 
reaction, 32 the condition of chemical reaction equilibrium, 



dG 

dx 



= 0, 



T.P 



defines the equilibrium fraction x e through 



ln^ + ln 



1 Xp 



RT 



(l-2xe) =0. 



(9) 



(10) 



If the excess Gibbs energy becomes large enough, phase sep- 
aration may occur, and a critical point may exist in the phase 
diagram. In the case of a regular solution, Eq. (8), the condi- 
tions for the critical point of liquid-liquid equilibrium, 



d 2 G \ 
dx 2 ) T p 



0. 



«9 3 G 



(11) 



T.P 



yield the critical composition x c = 1/2 and the condition 

T = HP) 

2R ' 



(12) 



If the line in the phase diagram given by Eq. (12) intersects the 
line given by the phase equilibrium condition lnK(T,P) = 0, 
a critical point exists at the intersection and a line of liquid- 
liquid phase equilibrium emanates from this point toward 
lower temperature. This is an energy-driven phase transition, 
like in the lattice-gas model. 
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FIG. 1. Density p, isobaric heat capacity Cp, isothermal compressibility Kt, and thermal expansivity dp computed from the mW model (points) compared 
with the results (curves) of the two-state approach for an athermal solution [Eq. (13)]. The isobar pressures are given in the density diagram; pressures and 
corresponding isobar colors are the same for the other plots. 



Athermal solution 



Clustering of water molecules 



In the case of an athermal solution, the enthalpy of mixing 
is zero and the non-ideality is associated with the excess en- 
tropy of mixing. With the simplest symmetric form of excess 
entropy, S E = —R(0x(l —x), we obtain 

G = G A +xG BA +RT[xlnx+(l -jc)ln(l -x) + cox(l -x)]. 

(13) 

The equilibrium fraction x e is given by 

ln^ + ln^^ + 0)(l -2x e ) = 0. (14) 
1 x e 

The critical-point conditions of Eq. (11) yield the critical value 



CO, 



(15) 



Thus, the critical pressure P c is the pressure at which the func- 
tion Co(P) reaches the value 2. The critical temperature T c fol- 
lows from the phase-equilibrium condition, In K(T C ,P C ) = 0. If 
the function co(P) remains below 2 for all pressures, a critical 
point does not exist and the mixture does not phase separate. 
If co is larger than 2 for a certain pressure, an entropy-driven 
phase transition exists. While real water shows a preference 
for the athermal two-state description, 15 the experimental data 
cannot exclude a contribution of energy in the non-ideal part 
of the Gibbs energy. 



Moore and Molinero 23 have demonstrated that four- 
coordinated molecules cluster in supercooled mW water, an 
idea originally proposed by Stanley and Teixeira. 33 The for- 
mation of clusters of molecules belonging to a single state 
reduces the number of configurations and decreases the mix- 
ing entropy. Clustering can be incorporated in the equation of 
state by dividing the ideal mixing entropy terms by the num- 
ber of molecules in a cluster N. For the athermal solution, 
Eq. (13), this yields 



G= G A +xG BA 
x 



RT 



Inx- 



1 -x 



ln(l —x) + cox(l —x) 



(16) 



Formally, this expression is identical to the free energy of a 
mixture of two polymers with (large) degrees of polymeriza- 
tion Ni and N2 in Flory-Huggins theory. 34 However, in our 
case these parameters are purely phenomenological. In prin- 
ciple, the cluster sizes Ni and N2 are functions of temperature 
and pressure. In this work, we use a constant N\ = N2 = N as 
a first approximation. Furthermore, one can imagine a mixed 
scenario in which the system combines both athermal-solution 
and regular-solution features. 
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FIG. 2. Fraction x of molecules in the low-density state, Solid curves: frac- 
tion x for the (a) two-state equation of state, Eq. (13); (b) two-state equation 
with account for hexamer clustering, Eq. (16) with N = 6. Dashed curves: 
fraction x obtained from simulations of mW water, calculated from the frac- 
tion of four-coordinated molecules f/s, as x = (f'4 — $)/(f\ — /4 1 ), to ac- 
count for fractions f\ and f? of four-coordinated molecules in the low- and 
high-temperature liquid, respectively. The inflection points on the curves are 
marked with circles (two-state equation) and squares (mW model). The data 
was collected by linearly quenching the temperature of the simulations at a 
rate of 10 K/ns. The data below the inflection point do not correspond to equi- 
librium states. 



Description of thermodynamic properties of the mW 
model 



Thermodynamic properties of the mW model, namely den- 
sity, isothermal compressibility, and heat capacity have been 
calculated from molecular dynamics simulations by Limmer 
and Chandler up to 228 MPa. 20 For this work, the expansion 
coefficient has also been calculated, and the calculations have 
been extended to higher pressures, as shown in Fig. 1. We 
apply the two-state thermodynamics of water to explain and 
reproduce the calculated properties and demonstrate that this 
model does not show liquid-liquid separation, at least in the 
range of pressures and temperatures studied. 

First, we have verified whether an "ideal-solution" two- 
state thermodynamics, for which the excess Gibbs energy is 
zero, can reproduce the thermodynamic properties of the mW 
model. Indeed, the competition between these two states, even 
without non-ideality, can cause the density maximum and the 
increase in the response functions. However, the thermody- 
namic properties are not the only data that should be matched. 
Moore and Molinero 23 have calculated the fraction of four- 
coordinated molecules fa, which is a fraction of molecules 
in a low-density state (Fig. 2). Specifically, fa is the fraction 
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FIG. 3. Pressure-temperature diagram. Circles: location of the computed 
property data of the mW model. Solid line: line at which InK = and the 
low-density fraction x — 1/2 for the two-state equation. Squares: location of 
inflection points of the low-density fraction x for the mW model (see Fig. 2). 
Dotted line: stability-limit temperature T s from the fit of the weak crystalliza- 
tion model to the mW data, Eq. (26). The dashed curve is a fit to the melting 
temperature T m of mW ice (uncertainty about ± 3 K), obtained from free 
energy calculations as described in Ref. 20. 



of molecules with four neighbors within the first coordination 
shell up to a certain cutoff radius r c . The exact value of fa de- 
pends on the value of r c , but the inflection point of the fa{T) 
curve is independent of r c and occurs at 7] = (201 ±2) K at 
atmospheric pressure. 23 The value of r c was 0.35 nm for all 
pressures in this study. This cutoff corresponds to the first sol- 
vation shell, as the position of the first minimum of the radial 
distribution function of mW water is 0.35 nm at 0.1 MPa and 
0.342 nm at 507 MPa. The fraction x of low-density liquid is 
related to fa by 



x(T) 



ft- 



ff 



(17) 



which accounts for the finite fraction ff of four-coordinated 
molecules in the high-temperature liquid, and the fraction 
ff < 1 of four-coordinated molecules in the low-temperature 
liquid. Both ff and ff are estimated by an extrapolation of 
the fraction fa to high and low temperature. Below 7], liquid 
mW water cannot be equilibrated without crystallization. 19 
Nevertheless, the fraction fa was also computed below 7] as in 
Ref. 23: in quenching simulations the temperature was varied 
linearly at 10 K/ns, the slowest rate that results in the vitrifica- 
tion of mW water. Liquid mW water can be equilibrated down 
to 7] at a cooling rate of 10 K/ns, but any property extracted 
from the quenching simulations at T < 7] may depend on the 
cooling rate. 

In the two-state thermodynamics, as given by Eq. (10) 
or Eq. (14), the inflection point of the fraction x at atmo- 
spheric pressure occurs where InK = 0, near the temperature 
Tq [Eq. (5)]. To match the low-density fraction of the mW 
model, To should be close to 7]. This is not the case for the 
ideal-solution two-state version, where To = (160 ±4) K, sig- 
nificantly below 7]. 

With a nonzero excess Gibbs energy, the two-state approach 
is able to reproduce both the thermodynamic properties and 
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the inflection point of the fraction in the mW model. As Fig. 3 
shows, the inflection points of the fraction of low-density liq- 
uid in mW simulations at different pressures form a straight 
line in the phase diagram. Since the two-state equation yields 
inflection points at the line InK = 0, the inflection points of 
the mW model can be matched by adopting suitable values 
of the slope a and intercept Tq of the InK — line [Eq. (4)]. 
When the InK = line is fixed in this way, the athermal- 
solution (entropy-driven non-ideality) version yields a bet- 
ter description than the regular-solution (energy-driven non- 
ideality) version; the sum of squared deviations of the regular- 
solution fit from the thermodynamic property data is about 
50% higher than that of the athermal-solution fit. More con- 
vincing evidence in favor of the athermal-solution approxi- 
mation comes from the direct computation of the enthalpy 
and entropy of mW water; see below. An approximation with 
a constant interaction parameter, w or ft), works reasonably 
well, and the description can be further improved by making 
the interaction parameter weakly dependent on pressure. In 
the case of the athermal-solution version, the two-state ther- 
modynamics for the mW model excludes liquid-liquid sepa- 
ration at any temperature or pressure, because the value of ft) 
is lower than 2 (for a pressure-independent ft), the optimum 
value is 1.61). In the case of the regular-solution version, the 
two-state thermodynamics for the mW model predicts liquid- 
liquid separation at a temperature below 147 K, which is out- 
side the kinetic limit of metastability of mW liquid water. 19 
Figure 1 shows the predicted thermodynamic properties for 
the athermal-solution case with a quadratic pressure depen- 
dence of ft), and Fig. 2a shows the low-density fraction. The 
numerical values of all parameters are given in the appendix. 
At 0.1 MPa, there is a systematic difference between com- 
pressibility values calculated from the two-state equation and 
those of the mW model. This difference can be decreased by 
including more background terms, but that would make the 
description more empirical. 

An improved description of the low-density fraction is ob- 
tained when clustering of water molecules is taken into ac- 
count in the equation of state, as in Eq. (16). When the number 
of molecules N ina cluster is taken as an adjustable parameter 
in the fit, the optimum value is 6.5, but the quality of the fit 
varies little for values of N between 4 and 10. For clusters of 
N = 6 molecules (hexamers), the fraction that results from the 
fit to the mW properties is shown in Fig. 2b. Because of the 
division of the ideal mixing entropy terms by N in Eq. (16), a 
smaller non-ideality term is sufficient. Indeed, for the fit with 
N = 6, the value of the interaction parameter ft) is 0.2, an order 
of magnitude smaller than in the case without cluster forma- 
tion. For such a small value of ft), the difference between a reg- 
ular and an athermal solution becomes unimportant, because 
the main contribution to the nonideality comes from clustering 
and thus is entropy driven. The work of Moore and Molinero 23 
shows that clustering of four-coordinated molecules increases 
with cooling, which suggests that the parameter Ni is tem- 
perature dependent. Such a temperature-dependent Ni would 
affect the values of the properties calculated with the two-state 
thermodynamics, particularly at low temperature. 

Analysis of the enthalpy of liquid water supports the con- 
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FIG. 4. Enthalpy AH (a) and entropy AS (b) of liquid mW water with respect 
to ice at 0.1 MPa (black curves, from Moore and Molinero") and their fits 
(dashed) according to Eqs. (18) and (19), respectively. Both AH and AS are 
computed at a cooling rate of 10 K/ns, which prevents crystallization, so that 
the values below 200 K do not correspond to an equilibrium state. The circle 
signals 7] = 201 K. These results support the modeling of mW water as an 
athermal mixture of two states. 



jecture that the excess free energy of mixing is almost entirely 
due to entropic effects. If water is represented as a "mixture," 
then its enthalpy can be written in terms of the partial mo- 
lar enthalpies of the low and high-density "components". Fig- 
ure 4 shows that the enthalpy of liquid water with respect to 
ice in mW water is very well represented by a sum of the 
weighted contributions of the two pure components, i.e., no 
excess enthalpy of mixing: 



AH 



[1.84x + 5.38(l-x)]kJ/mol. (18) 



Here we use the enthalpy with respect to ice, instead of the en- 
thalpy itself, to eliminate the trivial temperature dependence 
of the partial molar enthalpies of the components. This is 
possible in mW water because the heat capacity of the four- 
coordinated component is almost indistinguishable from that 
of mW ice, and the heat capacity of the high-temperature com- 
ponent is also quite close to that of ice, as a result of the 
monatomic nature of the model. Thus, the temperature depen- 
dence of the molar enthalpy difference between a component 
and ice is negligible. 

The enthalpy of the low-density component with respect 
to ice, 1.84 kJ/mol, is in good agreement with the (1.35 ± 
0.15) kJ/mol measured for LDA in experiments. 35 The en- 
thalpy of the high-density component with respect to ice is, 
unsurprisingly, the enthalpy of melting of ice (5.3 kJ/mol for 
mW, 6.0 kJ/mol in experiments). These results support the in- 
terpretation that there is no excess enthalpy of mixing contri- 
bution to the non-ideal excess free energy of liquid mW water. 

Interestingly, the entropy of liquid water with respect to ice 
can also be represented as a weighted sum of temperature- 
independent contributions from the two pure components 
(Fig. 4): 

AS = Si iq -S ice = [2.18jc+ 19.75(1 -*)] J/(Kmol), (19) 

where the entropy of the low-density component with respect 
to ice, 2.18 J/(Kmol), is in good agreement with the exper- 
imental value for LDA, (1.6 ± 1.0) J/(Kmol), 36 and the en- 
tropy of the high-density component with respect to ice is the 
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FIG. 5. Heat capacity of mW water in equilibrium (circles) and on hyper- 
quenching at 10 K/ns (black curve, computations by Moore and Molinero 19 ). 
The values below 200 K do not correspond to an equilibrium state. The dashed 
curve is the prediction of the two-state equation with hexamer clusters. 



entropy of melting (19.3 J/(Kmol) for mW water, 22 J/(Kmol) 
in experiment). The implication of Eq. (19) is that mW wa- 
ter is an "athermal solution," the excess (non-ideal) entropy 
of mixing is negative, and - remarkably - it essentially can- 
cels out the positive ideal entropy of mixing contribution. This 
cancellation strongly supports the idea of molecular cluster- 
ing and is similar to the thermodynamics of near-athermal 
mixtures of two high-molecular-weight polymers, where the 
entropy of mixing is very small. This property of mW water 
excludes the possibility of the regular-solution approximation. 

Figure 5 shows the isobaric heat capacity Cp calculated for 
liquid mW water down to a temperature below the crystalliza- 
tion temperature, which is made possible by hyperquenching 
at 10 K/ns. Below 200 K, the two-state thermodynamics pre- 
dicts a Cp that is lower than the value for hyperquenched mW 
water. This discrepancy may be caused by an oversimplicity 
of our equation of state, in particular the symmetric form of 
the excess entropy of mixing S = —Rcox(l —x) and the con- 
stant value of N. 

The two-state thermodynamics predicts maxima of the heat 
capacity and compressibility, and minima of the density and 
expansivity near the line InK = 0. These extrema are not ob- 
served in simulations of the equilibrium mW liquid, shown 
in Fig. 1, but they are seen when the system is driven out 
of equilibrium through fast cooling rates, as in Fig. 5. Fur- 
thermore, the two-state thermodynamics cannot predict the 
stability limit of the liquid phase or account for any pre- 
crystallization effects. 



III. WEAK CRYSTALLIZATION THEORY 

According to the theory of weak crystallization, 37-39 fluctua- 
tions of the translational order parameter cause corrections in 
the thermodynamic response functions close to the absolute 
stability limit of the liquid phase with respect to crystalliza- 
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FIG. 6. Fluctuation-renormalized distance to the stability-limit temperature 
A, given by Eq. (21), as a function of the mean-field distance to the stability 
limit Ao, for two values of j3 (solid curves). The dashed line corresponds to 
A = Aq and is shown as a reference. 



tion. The distance to the stability limit A is defined as 



A(T,P) 



T-T S {P) 



(20) 



where T is the temperature, P is the pressure, and T S (P) is 
the stability-limit temperature. According to this theory, the 
fluctuations of the translational (short-wavelength) order pa- 
rameter y/ renormalize the mean-field distance Ao = (T — 
r s MF )/r s MF between the temperature T and the mean-field ab- 
solute stability limit r s MF of the liquid phase: 



A = A + J3A 



-1/2 



(21) 



where A is the fluctuation-renormalized distance to the 
stability-limit temperature, and j3 is a molecular parameter, 
similar to the Ginzburg number that defines the validity of the 
mean-field approximation in the theory of critical phenomena. 
Solutions of Eq. (21) are shown in Fig. 6. Fluctuations of the 
translational order parameter stabilize the liquid phase, shift- 
ing the stability limit of the liquid below the mean-field value, 
meaning that Ao becomes negative at A tending to zero. This 
also means that the fluctuation part would not actually diverge. 
Thus the theory of weak crystallization requires J3A 1 72 <C Aq 
and Ao <C 1. As seen in Fig. 6, even the value of the coupling 
constant )3 =0.1 is already beyond the validity limit of the 
theory since it corresponds to the mean-field gap Ao = — 1. 

The contribution of order-parameter fluctuations to the iso- 
baric heat capacity Cp and isothermal compressibility Kj is 
proportional to A~ 3 / 2 , while the contribution to the density 
and entropy is ^A~ l l 2 . Accordingly, the fluctuation contribu- 
tion Sfi = j3 A 1 72 can be included in the chemical potential fl 



A=j3A 1 / 2 + A 1 '(f, J P), 



(22) 



where j5 ~ 2/3 (see the appendix) is the amplitude of the fluc- 
tuation part, and p. 1 is the regular (background) part of the 
chemical potential. The contributions from fluctuations are to 
be small compared to the background. The variables with a 
hat in Eq. (22) are dimensionless variables, defined as 



T = 



M = 



I 1 

RT, 



P = 



PVo 



(23) 



s() 




where T s o is the stability-limit temperature at zero pressure, 
R is the molar gas constant, and Vq — M x 10~ 3 m 3 /kg is an 
arbitrary reference constant for the molar volume, with M the 
molar mass of water. For our application of the weak crys- 
tallization theory to the mW model results, we represent the 
regular part fl 1 ' by a truncated Taylor-series expansion, 

(l'(t,P) = '£c mn t m P n , (24) 

similar to Eq. (7). We approximate the dependence of the 
stability-limit temperature on the pressure by a linear func- 
tion, 

f s (P) = l-a'P, (25) 

where the slope a' is an adjustable parameter. Expressions for 
the thermodynamic properties, found from derivatives of the 
chemical potential of Eq. (22), are given in the appendix. 

Fit to the mW data 

A fit of Eq. (22) to the mW data, shown in Fig. 7, results in a 
stability-limit temperature of 

T S (P)/K= 175 -0.093 (P/MPa), (26) 




0.0 I ■ ' ■ ' ■ ' J 

180 200 220 240 

Temperature (K) 

FIG. 8. Fluctuations of the crystallization order parameter from weak crystal- 
lization theory (i.e. short-wavelength density fluctuations) in arbitrary units 
as a function of temperature, fit with Eq. (28) (curve). 

which is shown in Fig. 3, and an amplitude of 40 

j3=2.3±0.4. (27) 

As shown in the appendix, the amplitude j3 is related to the 
coupling constant j5 as /3 ~ 2/3. To verify whether there is 
a microscopic manifestation of weak crystallization theory, 
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we have computed the normalized short-wavelength density 
fluctuations, corresponding to the fluctuations of the order pa- 
rameter \j/ from weak crystallization theory, at atmospheric 
pressure as a function of temperature. The quantity plotted in 
Fig. 8 is dimensionless but its magnitude is not defined unam- 
biguously. It is calculated by taking the expectation value over 
the Fourier transform of the density operator which results in 
the form (exp(igr) exp(— iqr)) where r is a particle position 
and q is the wavenumber. The wavenumber is chosen at the 
highest peak in the structure factor. The computation is per- 
formed for wN = 8000 particle system at a constant pressure 
of 0. 1 MPa, and each point is averaged over 10 ns. The fluctu- 
ations would diverge in the thermodynamic limit for a crystal, 
illustrating the broken symmetry. In the liquid phase, the or- 
der parameter (y) =0. Figure 8 shows the order parameter 
fluctuations together with a fit of the form 

(8 W ) 2 =b[(T-T s )/T s ]- l/ \ (28) 

where b is a constant proportional to the coupling constant j3 
in Eq. (21). From the fit we obtain b ~ 0.24 and T s ~ 181 K, 
close to the value of 175 K found above. Another way to es- 
timate T % is by extrapolating the surface tension between liq- 
uid and crystal as a function of supercooling, and finding the 
temperature at which it becomes zero. An extrapolation of the 
surface tension values of Limmer and Chandler 41 yields T s ~ 
170 K. The stability-limit temperature T s lies far below the 
temperature of maximum crystallization rate determined by 
Moore and Molinero, 19 which is 202 K at 0. 1 MPa and signals 
the lowest temperature at which liquid mW water can be equi- 
librated. This temperature difference should not be interpreted 
as a disagreement, because the temperature of maximum crys- 
tal growth is a kinetic quantity and does not correspond to a 
thermodynamic instability. 

While the weak crystallization model provides a reasonable 
description of the mW data as shown in Fig. 7, the ampli- 
tude /3 ~ 2j3 ~ 2 is not small compared to unity, contradicting 
Eq. (21), which is based on the assumption of small fluctua- 
tion corrections in the theory. The fit based on weak crystal- 
lization theory also implies that the density maximum at atmo- 
spheric pressure is caused by translational short-wavelength 
fluctuations of density. This is unphysical because such fluc- 
tuations should not have a significant effect far (100 K) from 
the stability limit of the liquid state, where the density maxi- 
mum is observed. Moreover, weak crystallization theory pre- 
dicts universal fluctuation-induced corrections to the regular 
behavior of thermodynamic properties and is equally applica- 
ble to all metastable fluids, not only to tetrahedral fluids that 
expand on freezing and exhibit a density maximum. 

That the weak-crystallization power laws provide a good 
empirical description of water's anomalies is not surprising. 
The properties of real supercooled water can also be described 
quite well with purely empirical power laws diverging along 
a certain line below the homogenous ice nucleation limit. 4 ' 42 
The fitted line of apparent singularities and the fitted values 
of the exponents are strongly correlated and cannot be ob- 
tained independently. The larger the exponent value, the fur- 
ther away from the homogenous nucleation this line is shifted. 



Adjustable amplitudes of the power laws would increase the 
ambiguity even further. 

However, we cannot ignore the fact that the fluctuations of 
the crystallization order parameter indeed increase upon su- 
percooling, as demonstrated in Fig. 8. While the relation be- 
tween the amplitude of this effect and the coupling constant j3 
of weak crystallization theory is unclear, it seems quite plausi- 
ble that the fluctuations also contribute to the thermodynamic 
anomalies. 



IV. CONCLUSIONS 

By fitting properties with equations of state based on differ- 
ent underlying assumptions, the consistency of such assump- 
tions with thermodynamics can be tested. In this way we have 
examined the origin of the anomalous behavior of the mW 
model in terms of a two-state model and a model based on 
weak crystallization. 

We have reproduced the anomalies of the mW model with 
a phenomenology based on a non-ideal mixture of two differ- 
ent molecular configurations in liquid water. While an ideal 
mixture of these states also generates enhancements in the re- 
sponse functions, agreement with the mW data on the low- 
density fraction unambiguously requires a non-ideal mixture. 
However, the nonideality is not strong enough to cause liquid- 
liquid phase separation in the mW model, in the wide range of 
pressures and temperatures of this study. An entropy-driven, 
athermal-solution-like nonideality gives a better description of 
the mW data than an energy-driven, regular-solution-like non- 
ideality. Incorporating the formation of clusters into the equa- 
tion of state results in a further improvement of the description 
of the low-density fraction in mW water. The thermodynamic 
treatment of the two-state model used here is independent of 
the details of the states considered. Microscopically, however, 
there is no unique projection for defining these states from the 
continuum distribution of configurations. 

A description which uses the power laws of weak crystal- 
lization theory succeeds in reproducing the calculated mW 
properties. The required-by-fit value of the coupling constant 
appears to be about 2, which contradicts the assumption of 
the theory where this constant is essentially a small param- 
eter. Finally, weak crystallization theory is based on the ex- 
istence of the liquid-state stability limit and accounts for the 
growing short-wave density fluctuations near that limit, which 
indeed are found in simulations. 19 ' 20 However, unlike the two- 
state approach, it cannot account for the data on the fraction 
of four-coordinated molecules. The need to consider fluctua- 
tion effects due to crystallization is consistent with previous 
work on the modulation of ice interfaces in real water. 41 This 
previous work demonstrated that such fluctuations are impor- 
tant in describing how crystallization is altered in confinement 
and suggests that such effects can be correctly accounted for 
within simple Gaussian corrections. Reconciling the anoma- 
lous thermodynamics that are well recovered by the two-state 
model, with the unavoidable effects of crystallization is an im- 
portant avenue to pursue in the future. 
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APPENDIX 

Expressions for thermodynamic properties in weak 
crystallization theory 

From derivatives of the chemical potential of Eq. (22), we ob- 
tain 



V = 

s = - 

C P = f 
Kr = - 



c)Pj T 2f 2 A 1 / 2 



D 

dfl 
df 

dl 
df 



pa'f 



dTjp 2f s AV2 ^' 



= T 



4 f 2 A 3/2 TT J , 

1 (dV\ 1 / pa' 2 f 2 _ Pa' 2 f _ 



1 fdV 



1 



pa'f 



pa' 



V \f) fjp V V 4f s 3 A 3 / 2 1 2f s 2A!/2 , fi fP 



(29) 
(30) 
(31) 

3 

(32) 

)■ 

(33) 



where V — V /Vb, S, Cp, kj and frp are the dimensionless vol- 
ume, entropy, isobaric heat capacity, isothermal compressibil- 
ity and the thermal expansivity, respectively. The subscripts of 
/i r indicate a partial derivative of fi 1 with respect to the sub- 
scripted quantities. 



Relation between the amplitude p and the coupling 
constant /3 in weak crystallization theory 

Consider the temperature-squared term in the dimensionless 
chemical potential ft, 



TABLE I. Parameters for the two-state equation of state, Eq. (13) 



Parameter 


Value 


Parameter 


Value 




T 

h) 




c 03 


1 .0 IZJ X 1U 


-4 


Po 


1UUU.U Kg m 


C05 


i cane w i n- 
— 1 .j lyj X 1U 


-7 


X 


2.6917 


C] 1 


8.4408 x 10~ 


-2 




1.6362 


C12 


-3.9163 x 10" 


-3 


CD, 


1.6777 


C13 


1.2385 x 10~ 


-4 


P 


2.3219 


C20 


- 1.7092 x 10° 




a 


0.039968 


C21 


-4.3879 x 10~ 


-2 


coi 


9.5440 x 10" 1 


C30 


4.0687 x 10- 


-1 


C02 


-5.0517 x 10~ 3 


C31 


6.1937 x 10- 


-2 




X 1 


= 0.82 






TABLE II. Parameters for the two-state equation, 


Eq. (16), withW = 


6 


Parameter 


Value 


Parameter 


Value 




To 


203.07 


<?03 


Z.O40 / X 1U 


-4 


Po 


1000.0 kg m" 3 


<?05 


— j. Lojo X 1U 


-7 


X 


1.529 


cu 


3.9710 x 10- 


-2 


COo 


0.19523 


C\2 


-6.4543 x 10- 


-4 


CD, 


0.19856 


C13 


4.6055 x 10- 


-5 


Pi 


1.1637 


C20 


-1.8630 x 10° 




a 


0.039968 


C21 


-3.2135 x 10- 


-2 


CQ\ 


9.8530 x 10" 1 


C30 


3.3420 x 10- 


-1 


C02 


-8.1978 x 10~ 3 


C31 


5.7351 x 10- 


-2 




X 2 


= 0.72 







Considering the A term, using Eq. (21) gives 



2j3A A 



1/2. 



J3 2 A 



= A 2 , + 2j3A 1 / 2 - i S 2 A 



2a-1 



(35) 



The last term can be neglected because according to the the- 
ory ft <C 1 . Since this is not the case in our attempt to apply 
weak crystallization theory to the mW data, such a description 
becomes purely empirical. 

The first term on the right-hand side Aq is part of the regu- 
lar part of the chemical potential, and the second term 2/3A 1 / 2 
leads to the fluctuation contribution dp. in the chemical poten- 
tial of 



<5£ = 2cf s 2 pA 



1/2 



(36) 



Comparing this with Eq. (22) gives the relation between j3 and 
P. 



P=2cf 2 p. 



(37) 



cf 2 = cf 2 (A+ l) 2 = cf 2 (A 2 + 2A+l) 



(34) 



For our fit to the mW data, c ~ 1 and T s ~ 1, so Eq. (37) gives 
P-2p. (38) 

Parameter values 

Parameters for the two-state equation of state, Eq. (13), are 
listed in Table I, and those for the fit of the equation with clus- 
ters of six molecules, Eq. (16), are given in Table II. The value 



10 



TABLE III. Parameters for the fit of weak crystallization theory, Eq. (22) 



Parameter 


Value 




Parameter 


Value 




7k) 


175.28 K 




Cll 


4.2389 x 10" 


-I 


M/V 


1000.0 kg rrr 3 




C\2 


-3.3672 x 10" 


-4 


jS 


2.3254 




C13 


-1.3649 x 10" 


-4 


d 


0.042825 




C20 


3.4178 x 10" 


-1 


coi 


6.1915 x 10" 


-l 


C21 


-2.3514 x 10" 


-1 


C02 


-8.2724 x 10" 


-3 


C30 


-2.2560 x 10" 


-1 


C03 


3.6026 x 10" 


-4 


C31 


4.9163 x 10" 


-2 


c 05 


-6.6486 x 10" 


-7 









t = 0.72 

of the interaction parameter (O depends quadratically on the 
pressure, 



co{P) = (o x - {(o x - coo) [(P- A)/A] : 



(39) 



where ©o is the value of (0 at P = 0, C0i is the maximum value 
of o, and Pi is the dimensionless pressure at which the max- 
imum occurs. The parameter values for the fit of weak crys- 
tallization theory to the mW data are given in Table III. All 
tables give % 2 , the reduced sum of squared deviations of the 
fit from the data. 
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